clearvars -except b
load datafromstata tprob_m v_crime Lv_crime  tprob_t LLv_crime

tempind = (tprob_t>2);
tprob_m = tprob_m(tempind);
v_crime = v_crime(tempind);
Lv_crime = Lv_crime(tempind);
tprob_t = tprob_t(tempind);
LLv_crime = LLv_crime(tempind);


num_coef = 4;

M = max(tprob_m);
pxest = zeros(M,num_coef);sigrhoest=zeros(M,1);

for m = 1:M
    index = find(tprob_m==m);
    cons = ones(length(index),1);
    [~ ,~,r,~,~] = regress(v_crime(index), [cons Lv_crime(index) tprob_t(index) LLv_crime(index)]);
    temp = randsample(r,length(r),'true');
    [temp ,~,r,~,stats] = regress(v_crime(index)+temp, [cons Lv_crime(index) tprob_t(index) LLv_crime(index)]);
    pxest(m,:) = temp;
end
pxest = pxest';


T = 7; 

yvals = [1:1:19]';
y=1;
gamma = ones(M,1);
beta = 0.95;
xl = 2;

v_crime_bar_big = 999*ones(M,2);
for val = 1:2
    x = val-1;
    v_crime_bar = x*ones(M,1);
extp=x*ones(M,1);
extpL=xl*ones(M,1);
for m=1:M
B = 1;
for t=2:7 
extpL(m) =  extp(m);    
extp(m) =  pxest(1,m) + pxest(2,m)*extp(m)  + pxest(3,m)*(yvals(y)+t-1) + pxest(4,m)*extpL(m) ;


temp = beta^(t-1);
temp2 = (beta^(t-1))*( pxest(2,m)^(t-1));
temp3 = (beta^(t-1))*( (pxest(2,m)^(t-2))*(t-1));
B = B+temp;
v_crime_bar(m) = v_crime_bar(m) + temp*extp(m);
gamma(m) = gamma(m) + temp2;
end
end
v_crime_bar = v_crime_bar/B;
v_crime_bar_big(:,val) = v_crime_bar;
end

gammaar2 = v_crime_bar_big(:,2) - v_crime_bar_big(:,1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load est_res_s2 rentd  M con_m
load datafromstata gamma 

statmwtp = rentd;

picmwtpar2 = 0*rentd;
for m = 1:M
    index = find(con_m==m);
    picmwtpar2(index) = statmwtp(index)/gammaar2(m);
end

for m = 1:M
    index = find(con_m==m);
    res2(m,:) = [ mean(picmwtpar2(index))];
end

rob_AR2 = res2;
save rob_AR2 rob_AR2


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clearvars -except b

load tprob_data 


M = max(tprob_m);
pxest2 = 99999+zeros(M,max(tprob_t));


yvals = [1:1:19]';

for m = 1:M
    for y=1:length(yvals);
    index = find(tprob_m==m & tprob_t <= y);
    if ~isempty(index)
    cons = ones(length(index),1);
    [~ ,~,r,~,~] = regress(v_crime(index), [cons Lv_crime(index) tprob_t(index)]);
    temp = randsample(r,length(r),'true');
    [temp ,~,r,~,stats] = regress(v_crime(index)+temp, [cons Lv_crime(index) tprob_t(index)]);
    pxest2(m,y) = temp(2);
    end
    end
end



T = 7; beta = .95;

gammaae = zeros(M,19);
for m = 1:M
    %m
    for y=1:length(yvals);
    if pxest2(m,y) < 99999
    gamma = 1;
    B = 1;
    for t=2:T
    temp = beta^(t-1);
    temp2 = (beta^(t-1))*( pxest2(m,y)^(t-1));
    B = B+temp;
    gamma = gamma + temp2;
    end
    gamma = gamma/B;
    gammaae(m,y) = gamma;          
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load est_res_s2 rentd  M con_m t
load datafromstata gamma 

statmwtp = rentd;

picmwtpae = 0*rentd;  
for m = 1:M
    for y=3:length(yvals);
        index = find(con_m==m & y==t);
        picmwtpae(index) = statmwtp(index)/gammaae(m,y);
    end  
end

for m = 1:M
    index2 = find(con_m==m & t >= 3);
    res2(m,:) =  mean(picmwtpae(index2));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rob_re = res2;
save rob_re rob_re

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clearvars -except b
load tprob_data 



num_coef = 4;

M = max(tprob_m);
pxest = zeros(M,num_coef);sigrhoest=zeros(M,1);

for m = 1:M
    index = find(tprob_m==m);
    cons = ones(length(index),1);
    [~ ,~,r,~,stats] = regress(v_crime(index), [cons Lv_crime(index) tprob_t(index) Lpm10(index)]);
    temp = randsample(r,length(r),'true');
    [temp ,~,r,~,stats] = regress(v_crime(index)+temp, [cons Lv_crime(index) tprob_t(index) Lpm10(index)]);
    pxest(m,:) = temp;
end
pxest = pxest';

T= 7; beta = 0.95;

gamma_tprob =  gammafunc(T,beta, M, pxest);

load est_res_s2 rentd  M con_m 
load datafromstata gamma

statmwtp = rentd;

picmwtpr1 = 0*rentd; 
for m = 1:M
    index = find(con_m==m);
    picmwtpr1(index) = statmwtp(index)/gamma_tprob(m,1);
end

for m = 1:M
    index = find(con_m==m);
    res2(m,:) = mean(picmwtpr1(index));
end

rob_tprob = res2;
save rob_tprob rob_tprob 
   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
robust_T_beta 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clearvars -except b
load tprob_data 

poscon = max(max([v_crime Lv_crime]));
v_crime = log(poscon*1.1 - v_crime);
Lv_crime = log(poscon*1.1 - Lv_crime);


num_coef = 3;

M = max(tprob_m);
pxest = zeros(M,num_coef);sigrhoest=zeros(M,1);


for m = 1:M
    index = find(tprob_m==m);
    cons = ones(length(index),1);
    [~ ,~,r,~,~] = regress(v_crime(index), [cons Lv_crime(index) tprob_t(index) ]);
    temp = randsample(r,length(r),'true');
    [temp ,~,r,~,stats] = regress(v_crime(index)+temp, [cons Lv_crime(index) tprob_t(index) ]);
    pxest(m,:) = temp;
end
pxest = pxest';

T = 7; 

yvals = [1:1:19]';
y=1;
gamma = ones(M,1);
beta = 0.95;
xl = 2;


v_crime_bar_big = 999*ones(M,2);
for val = 1:2
    x = val-1;
    v_crime_bar = x*ones(M,1);
extp=x*ones(M,1);
extpL=xl*ones(M,1);
for m=1:M
B = 1;
for t=2:7 
extpL(m) =  extp(m);    
extp(m) =  pxest(1,m) + pxest(2,m)*extp(m)  + pxest(3,m)*(yvals(y)+t-1) ;

temp = beta^(t-1);
temp2 = (beta^(t-1))*( pxest(2,m)^(t-1));
temp3 = (beta^(t-1))*( (pxest(2,m)^(t-2))*(t-1));
B = B+temp;
v_crime_bar(m) = v_crime_bar(m) + temp*extp(m);
gamma(m) = gamma(m) + temp2;
end
end
v_crime_bar = v_crime_bar/B;
v_crime_bar_big(:,val) = v_crime_bar;
end

gammalog = v_crime_bar_big(:,2) - v_crime_bar_big(:,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load est_res_s2 rentd  M con_m 
load datafromstata gamma 

statmwtp = rentd;
picmwtplog = 0*rentd;
for m = 1:M
    index = find(con_m==m);
    picmwtplog(index) = statmwtp(index)/gammalog(m);
end

for m = 1:M
    index = find(con_m==m);
    res2(m,:) = [ mean(picmwtplog(index))];
end

rob_log = res2;
save rob_log rob_log



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clearvars -except b
load rob_AR2
load rob_re 
load rob_tprob
load rob_Tb
load rob_log
rob_res = [rob_base  rob_AR2 rob_re rob_tprob rob_Tb rob_log]';

save rob_res rob_res